Quadratwurzel ziehen mittels Assembler

Im Rahmen des Ray Casting bin ich auf eine weitere Hürde gestossen: Zur Berechnung der Distanz zwischen Player und Wand/Hürde brauche ich die Quadratwurzel. Mehrere Möglichkeiten fallen mir nach langer Überlegung dazu ein: Tabellen, bei einem großen Playfield wird das aber recht speicherintensiv... oder berechnen! Die mir bekannte schriftliche Berechnung ist meines Erachtens schwierig umzusetzen, daher habe ich weiter recherchiert und bin auf eine andere einfachere Lösung gestossen.... es geht um die Heron-Methode!

Als Beispiel nehme ich mal die Wurzel aus 12 - laut Taschenrechner = 3,464 Man definiert ein Rechteck mit 1x12 und nähert sich dann dem Quadrat an. Zunächst wird ein x0 gebildet: x0 = (12+1)/2 = 6,5 Dann x1: x1 = (6,5+12/6,5) /2 = 4,173 Dann x2: x2 = (4,173 + 12/4,173) / 2 = 3,52 Dann x3: x3 = (3,52 + 12/3,52) / 2 = 3,465

Mit jedem Iterationsschritt wird das Ergebnis genauer.

Wenn mir das Ergebnis etwas gröber reicht und ich die Kommastellen weglasse kommt folgendes raus:

x0 = 6,5 = 6 x1 = 4 x2 = 3,5 = 3 x3 = 3,5 = 3

Edit: Bei hohen Zahlen sollte man eine Iteration mehr nehmen!! So liefert die Methode bei Wurzel aus 134 (=11,5758) nach 3 Schritten mit Abrundung das Ergebnis 13 und nach einem Schritt mehr kommt man auf 11,5 also abgerundet auf 11 und aufgerundet auf 12!

Als nächstes gilt es diesen Heronalgorithmus in Assembler umzusetzen....

Gibt es noch einfachere Möglichkeiten oder Umsetzungen????

CORDIC kannte ich noch nicht - sehr interessant / vielen Dank - scheint aber eher langsam zu sein und setzt doch ein wenig mathematisches Verständnis voraus (Vektorrotation etc.) - das dann noch in Assembler umzusetzen ist sicherlich sportlich.

Gute Idee - da werde ich ihn unbedingt mal drauf ansprechen und ne kleine Mathestunde von ihm erbeten .. bin ja selber so eine Art Hobbymathematiker 😉

sehr cool ... welchem Algorithmus liegt dem zugrunde??

Hallo Lars, das geht natürlich im Prinzip auch - und ich denke bei Pacman ist das super schnell und sinnvoll - die Näherung (x1-x2) + (y1-y2)/2 wenn y1-y2 die kurze Seite ist, ist deutlich genauer. Ich denke beim Raytracing sollte es etwas genauer sein, da vielleicht sonst eine Verzerrung auftritt.... ich weiß es aber noch nicht. Ich werde beides testen und dann posten! LG Peter

Irgendwer schrieb: Meine Routine aus "Disc o' Pop":

``` ; UWORD fastcall ; SM_Sqrt(ULONG lValue_); ; .define result sreg .define input ptr4 .define T ptr1

.export _SM_Sqrt .proc _SM_Sqrt

  sta input
  stx input+1
  lda sreg
  sta input+2
  lda sreg+1
  sta input+3

  lda #0
    sta result        ; R=0
    sta result+1
    sta input+4
    ;sta T+15        ; (T+1) is zero until last iteration; (T+0) is always 0

    clc
    ldy #14          ; 15 iterations (14-->0) + last iteration

loopsq:
lda result ; (2*R+D) LSR 1; actually: R+(D LSR 1) ora stablo,y ; using ORA instead of ADC is ok because the bit to be set sta T+2 ; will have not been affected yet lda result+1 ora stabhi,y sta T+3 bcs skip0 ; takes care of large numbers; if set, input>T

    lda input+3
    cmp T+3
    bcc skip1        ; T <= input    (branch if T>input)
    bne skip0
    lda input+2
    cmp T+2
    bcc skip1

skip0: ;sec lda input+2 ; input=input-T sbc T+2 sta input+2 lda input+3 sbc T+3 sta input+3 lda result ; R=R+D ora stablo+1,y sta result lda result+1 ora stabhi+1,y sta result+1 skip1: asl input ; input=input2 rol input+1 rol input+2 rol input+3 dey ; implicit: D=D/2, by the decrement of .Y bpl loopsq lastiter: ; code for last iteration bcs skp0 ; takes care of large numbers; if set, input>T ; during last iteration D=1, so [(2R+D) LSR 1] makes D the MSB of T+1 lda input+3 cmp result+1 ; (T+3) = resultHI bcc skp1 ; T <= input branch if T>input bne skp0 lda input+2 cmp result ; (T+2) = resultLO bcc skp1 bne skp0 lda input+1 cmp #$80 ; value of (T+1) during last iteration bcc skp1 skp0: ;sec lda input+1 sbc #$80 ; (T+1) during last iteration sta input+1 lda input+2 ; input=input-T sbc result ; (T+2) sta input+2 lda input+3 sbc result+1 ; (T+3) sta input+3 inc result ; R=R+D with D=1 skp1:
lda result ldx result+1 rts

.endproc

.data stabhi: .byte 0,0,0,0,0,0,0,0 stablo: .byte $01,$02,$04,$08,$10,$20,$40,$80 .byte 0,0,0,0,0,0,0,0 ```

welcher Idee bzw. Routine/Algorithmus steckt denn dahinter???

... ja Heron ... leider englisch und nicht super kommentiert...

Ich werde hier mal für MADS eine Heron-Routine machen, die wir dann auch ins Wiki stellen können...

.... ich habe das Wurzel ziehen mittels Heron-Iteration zunächst mal in Basic umgesetzt und geschaut, wie genau das Ergebnis nach x-Iterationsschritten aussieht. Wenn man im 8-Bit Bereich bleibt, dann reichen 6-7 Schritte völlig aus für ein sehr genaues Ergebnis. Geht man in den 16-Bit Bereich und will z.B. von 65000 die Quadratwurzel ziehen, dann braucht man 10 Durchläufe.

Hier zunächst der Basic-Code:

5   Graphics 0
10  Input A
20  B=1
30  C=A
40  For F=1 To 6
50  C=C+B
60  C=C/2
70  B=A/C
80  Next F
90  Print C

Mit B=1 und C=A definiere ich das erste Rechteck C x B (9x1) - nach Heron wird nun Schritt für Schritt die Seite C verkleinert und B vergrößert, bis beide Seiten sich annähern. Wenn beide Seiten annähernd gleich sind, dann steht das Ergebnis fest - die eine Seite (hier C) hat dann das Ergebnis Wurzel von A! Schon schlau die alten Griechen....

Bevor ich nun den Code in Assembler umsetzen, brauche ich noch eine Routine fürs Teilen!

... und hier jetzt der Assemblercode für die Quadratwurzel im Wertebereich 0-255 (8 Bit):

Wurzel aus var1 -> Ergebnis in var3

             .proc qwurzel  ;von var1
             ldx #0
             ldy #6
             lda var1
             sta var3
             lda #1
             sta var2
    lp1
            clc
            lda var3
            adc var2
            sta var3
            lsr var3
            lda var1
            ;------------ hier Divisionsroutine
            ldx #255
            sec
    lp2     
            sbc var3
            inx
            bcs lp2
            stx var2
            ;----------- Division Ende
            dey
            bpl lp1
            rts
            .endp

y muss auf Null sein - stimmt y rückwärts macht Sinn und ist schneller! Habe es gerade oben geändert!! Tak!!

als nächster Schritt wäre jetzt noch interessant, die Quadratwurzel auf eine 16 Bit- Zahl zu erweitern, dann bebommt man ruck zuck auf noch Pythagoras... Freiwillige vor 😉

tolle Routine! Ich finde es in der Mathematik und in der Informatik immer wieder faszinierend, wie viele verschiedene Ansätze zur gleichen Lösung führen!!

Ich habe jetzt die Routine für die Quadratwurzel einer 16-Bit Zahl nochmal analysiert und getestet und entsprechend der Variablen angepasst. Diese Routine kommt auch ins Atari-Wiki und ist hinsichtlich Schnelligkeit und Kreativität ein echter Schatz!!! In Gedenken an Lee Davison, der leider schon früh gestorben ist....

In der Routine wird die Quadratwurzel von var1/var1+1 gebildet. Das Ergebnis ist dann in var2, der Rest in var3. Ich habe versucht es zu durchdringen - es ist mir nicht vollständig gelungen! Insgesamt wurde hier ein schriftliches Wurzelziehen auf 2-Bit Ebene umgesetzt (daher die vielen rol und asl Anweisungen). Die Routine erscheint wirklich stark optimiert und ist entsprechend auch schön kompakt und schnell! Ich denke Lee hat da lange dran getüftelt und hat meinen uneingeschränkten Respekt dafür!!

        .proc qwurzel16bit
        ; Quadratwurzel 16 Bit
        ; nach Lee Davison
        ; Quadratwurzel var1/var1+1

        lda #0
        sta var2
        sta var3
        ldx #8
lp1
        sec
        lda var1+1
        sbc #64
        tay
        lda var3
        sbc var2
        bcc lp2
        sty var1+1
        sta var3
lp2
        rol var2
        asl var1
        rol var1+1
        rol var3
        asl var1
        rol var1+1
        rol var3
        dex
        bne lp1
        rts
        .endp

.... vielleicht ist hier auch der CORDIC Algorithmus verbaut .... mit diesem werde ich mich auch noch demnächst etwas näher beschäftigen!